#include <boost/filesystem.hpp>
#include <boost/mpi.hpp>
#include <cmath>
#include <eigen3/Eigen/Dense>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
//#include "checkTwoPi.h"
//#include "checkVsTok.h"
// #include "../SRC/DF_Common/const_nucl.h"
// #include "../SRC/DF_Common/coupling_coeff.h"
// #include "../SRC/DF_Common/file_manage.h"
// #include "../SRC/DF_Common/get_input_file.h"
// #include "../SRC/DF_Common/package.h"
// #include "../SRC/DF_Common/serial_control.h"
// #include "../SRC/DF_Common/time_count.h"
// #include "../SRC/DF_Common/basis_state.h"
// #include "../SRC/DF_Mpi_Sep/antisym_coeff.h"
// #include "../SRC/DF_Mpi_Sep/cont_term.h"
// #include "../SRC/DF_Mpi_Sep/info_wave.h"
// #include "../SRC/DF_Mpi_Sep/ope_term.h"
// #include "../SRC/DF_Mpi_Sep/tpe_term.h"
// #include "../SRC/DF_Mpi_Sep/trans_coeff.h"
// #include "../SRC/DF_Mpi_Sep/wavefun_coeff.h"
#include "../SRC/DF_Mpi_Sep/main_control.h"

// using std::cin;
using std::cin;
using std::cout;
using std::endl;
using std::string;
using std::vector;
namespace bf = boost::filesystem;
// using namespace std;

int main() {

  const std::string path_serial = "serial_dat";
  // bf::path path(path_s);
  if (!bf::exists(path_serial)) {
    bf::create_directory(path_serial);
  }
  const std::string path_out = "serial_dat/Output";
  // bf::path path(path_s);
  if (!bf::exists(path_out)) {
    bf::create_directory(path_out);
  }

  mpi::environment env_main;
  mpi::communicator MainCom;
  Eigen::initParallel();
  vector<int> mpi_run_flag_gather(MainCom.size());
  int mpi_run_flag = 1;
  // omp_set_num_threads(15);

  LEC_CONST LECs;
  if (MainCom.rank() == 0) {
    cout << " ---- LECs : ----" << endl;
    cout << "\t c1 = " << LECs.LEC_C1 << endl;
    cout << "\t c3 = " << LECs.LEC_C3 << endl;
    cout << "\t c4 = " << LECs.LEC_C4 << endl;
    cout << "\t cD = " << LECs.LEC_CD << endl;
    cout << "\t cE = " << LECs.LEC_CE << endl;
    cout << " ----------------" << endl;
  }
  Time_Count TC_main;
  TC_main.start();

  Get_Info get_info;

  string Input_src = "data/Input/DF.Input";

  get_info.init(Input_src);
  get_info.build();

  // exit(0);

  JaJbJ JaJbJ_t;
  int j2a_max = 9;
  int j2a_max_cg = 25;
  int J2max, J2min, J2num;
  Serial_SL ser_t;
  File_Manage File_M;
  vector<int> s_flag(9);
  char u1[256]; // unuse
  if (MainCom.rank() == 0) {
    std::ifstream ser_file;
    ser_file.open("data/Input/Step_Control.Input");

    // ser_file.getline(u1, 256, '\n');
    for (int i = 0; i < 9; i++) {

      int A = -1;
      ser_file >> A;
      ser_file.getline(u1, 256, '\n');
      // cout << A << "\t** " << endl;
      s_flag[i] = A;
    }
    ser_file.close();
  }
  boost::mpi::broadcast(MainCom, s_flag, 0);

  int serial_flag_SC_cal;   // 1: do the calculation other: load
  int serial_flag_SC_save;  //
  int serial_flag_WT_cal;   // 1: do the calculation other: load
  int serial_flag_WT_save;  //
  int serial_flag_OCP_cal;  // 1: do the calculation other: load
  int serial_flag_OCP_save; //
  int State_Coeff_Break;
  int Wave_Tpe_Break;
  int Test_Press;
  serial_flag_SC_cal = s_flag[0];   // 1:save other:load
  serial_flag_SC_save = s_flag[1];  //
  serial_flag_WT_cal = s_flag[2];   // 1:save other:load
  serial_flag_WT_save = s_flag[3];  //
  serial_flag_OCP_cal = s_flag[4];  // 1:save other:load
  serial_flag_OCP_save = s_flag[5]; //
  State_Coeff_Break = s_flag[6];
  Wave_Tpe_Break = s_flag[7];
  Test_Press = s_flag[8];
  if (MainCom.rank() == 0) {
    cout << "serial_flag_SC_cal (1:yes)   : " << serial_flag_SC_cal << endl;
    cout << "serial_flag_SC_save (1:yes)  : " << serial_flag_SC_save << endl;
    cout << "serial_flag_WT_cal (1:yes)   : " << serial_flag_WT_cal << endl;
    cout << "serial_flag_WT_save (1:yes)  : " << serial_flag_WT_save << endl;
    cout << "serial_flag_OCP_cal (1:yes)  : " << serial_flag_OCP_cal << endl;
    cout << "serial_flag_OCP_save (1:yes) : " << serial_flag_OCP_save << endl;
    cout << " --- ---- ---" << endl;
    cout << "State_Coeff_Break (1:break) : " << State_Coeff_Break << endl;
    cout << "Wave_Tpe_Break (1:break)    : " << Wave_Tpe_Break << endl;
    cout << "Test press     (1:test)    : " << Test_Press << endl;
  }

  if (serial_flag_SC_save == 1 && MainCom.rank() == 0) {
    const std::string path_s = "serial_dat/Output";
    // bf::path path(path_s);
    if (!bf::exists(path_s)) {
      bf::create_directory(path_s);
    }
    string f_dst = "serial_dat/Output/DF.Input";
    File_M.copy(Input_src, f_dst);
  }

  CG_coeff cg_x;
  SixJ_coeff sixJ_x;
  Binom bin_t;

  if (serial_flag_SC_cal == 1) {
    if (MainCom.rank() == 0)
      cout << "\t ---- ---- Start Part : Wave & Tpe ---- ---- " << endl;
    cg_x.init(10);

    sixJ_x.init(12);

    int n_max = 100;
    bin_t.init(n_max);
    if (serial_flag_SC_save == 1 && MainCom.rank() == 0) {
      const std::string path_s = "serial_dat/Output";
      // bf::path path(path_s);
      if (!bf::exists(path_s)) {
        bf::create_directory(path_s);
      }
      ser_t.save(cg_x, "serial_dat/Output/cg_serial.dat");
      ser_t.save(sixJ_x, "serial_dat/Output/sixJ_serial.dat");
      ser_t.save(bin_t, "serial_dat/Output/binom_serial.dat");
    }
  } else {
    ser_t.load(cg_x, "serial_dat/Output/cg_serial.dat");
    ser_t.load(sixJ_x, "serial_dat/Output/sixJ_serial.dat");
    ser_t.load(bin_t, "serial_dat/Output/binom_serial.dat");
    if (MainCom.rank() == 0)
      cout << "Finished Load cg_x, sixJ_x & bin_t" << endl;
  }
  // exit(0);
  NinJ_coeff ninJ_x(sixJ_x);
  if (serial_flag_SC_cal == 1) {
    ninJ_x.init(5);
    if (serial_flag_SC_save == 1 && MainCom.rank() == 0) {
      ser_t.save(ninJ_x, "serial_dat/Output/ninJ_serial.dat");
    }
  } else {
    ser_t.load(ninJ_x, "serial_dat/Output/ninJ_serial.dat");
    if (MainCom.rank() == 0)
      cout << "Finished Load ninJ_x" << endl;
  }
  // exit(0);
  Hobra hobra_t(bin_t, cg_x, ninJ_x);

  SingleP single_t;
  if (serial_flag_SC_cal == 1) {
    // const std::string path_sp = "data/Input/SP.Input";
    // const std::string path_sp_E = "data/Input/SP_E.Input";
    // bf::path path(path_s);
    // if (bf::exists(path_sp)) {
    //   single_t.build_T(path_sp);
    //   cout << "Build SP using data/Input/SP.Input " << endl;
    // } else if (bf::exists(path_sp_E)) {
    //   single_t.build_T_E(path_sp_E);
    //   cout << "Build SP using data/Input/SP_E.Input " << endl;
    // }
    if (get_info.sp_info.Type == "Read") {
      single_t.build_T(get_info.sp_info.sp_vec);
    } else if (get_info.sp_info.Type == "E_max") {
      single_t.build_T_E(get_info.sp_info.E_max);
    } else {
      cerr << " Wrong @ build SP " << endl;
      cerr << "\t get_info.sp_info.Type : " << get_info.sp_info.Type << endl;
      exit(0);
    }

    single_t.build_Tz();
    single_t.print_Tz();
    // exit(0);
    // single_t.build_Tz(0,1,1,1);
    // single_t.build_TzM(0,1,1,1);
    if (MainCom.rank() == 0 && serial_flag_SC_save == 1) {
      single_t.print("serial_dat/Output/SP.dat");
      single_t.print_Tz("serial_dat/Output/SP_Tz.dat");

      // single_t.print_TzM("data/SP_TzM.dat");

      ser_t.save(single_t, "serial_dat/Output/SP_serial.dat");
    }
  } else {
    ser_t.load(single_t, "serial_dat/Output/SP_serial.dat");
    if (MainCom.rank() == 0)
      cout << "Finished Load single_t" << endl;
  }

  // ------   ------ Start to build TP States ------ ------ //
  TripleP triple_t;
  vector<vector<int>> E_Jac_conf_t;

  vector<int> abc_vec(6);
  int TP_allow_file_num;
  int TP_E_max;
  // if (MainCom.rank() == 0) {
  //   std::ifstream TPinput_file;
  //   TPinput_file.open("data/Input/TP.Input");
  //
  //   for (int i = 0; i < 6; i++) {
  //
  //     int A = -1;
  //     TPinput_file >> A;
  //     // cout << A << "\t** " << endl;
  //     abc_vec[i] = A;
  //     TPinput_file.getline(u1, 256, '\n');
  //   }
  //   TPinput_file >> TP_allow_file_num;
  //   TPinput_file >> TP_E_max;
  //   TPinput_file.close();
  //   string f_src = "data/Input/TP.Input";
  //   string f_dst = "serial_dat/Output/TP.Input";
  //   File_M.copy(f_src, f_dst);
  // }
  // int master_node = 0;
  // boost::mpi::broadcast(MainCom, abc_vec, master_node);
  // boost::mpi::broadcast(MainCom, TP_allow_file_num, master_node);
  // boost::mpi::broadcast(MainCom, TP_E_max, master_node);
  //

  vector<vector<vector<int>>> TP_allow_gather;
  vector<int> E_cut_gather;
  vector<int> Same_Particles_gather;
  vector<vector<int>> E_sign_gather;

  TP_allow_file_num = get_info.tp_al_info.size();
  for (int i = 0; i < TP_allow_file_num; i++) {
    vector<vector<int>> TP_allow_vec(2);
    for (int j = 0; j < 2; j++) {
      TP_allow_vec[j].resize(6);
      for (int k = 0; k < 6; k++) {
        TP_allow_vec[j][k] = get_info.tp_al_info[i].tp_al_vec[k][j];
      }
    }
    vector<int> E_sign(6);
    for (int k = 0; k < 6; k++) {
      E_sign[k] = get_info.tp_al_info[i].tp_al_vec[k][2];
    }
    int Same_Particles_Flag = get_info.tp_al_info[i].Same_particle_Flag;
    int E_cut = get_info.tp_al_info[i].E_cut;
    TP_allow_gather.push_back(TP_allow_vec);
    E_cut_gather.push_back(E_cut);
    E_sign_gather.push_back(E_sign);
    Same_Particles_gather.push_back(Same_Particles_Flag);
    single_t.print(TP_allow_vec);
  }
  // exit(0);

  //
  // int TP_allow_file_num = 2;
  // if (MainCom.rank() == 0) {
  //   std::ifstream TP_Allow_file;
  //   for (int i = 0; i < TP_allow_file_num; i++) {
  //     vector<vector<int>> TP_allow_vec;
  //     TP_allow_vec.resize(2);
  //     for (int j = 0; j < 2; j++) {
  //       TP_allow_vec[j].resize(6);
  //     }
  //     vector<int> E_sign(6);
  //     int Same_Particles_Flag;
  //     int E_cut;
  //     int ii = i;
  //     string filename = "data/Input/TP_allow_" + std::to_string(ii);
  //     filename += ".Input";
  //     string filename_copy = "serial_dat/Output/TP_allow_" +
  //     std::to_string(ii); filename_copy += ".Input"; File_M.copy(filename,
  //     filename_copy); TP_Allow_file.open(filename); TP_Allow_file >> E_cut;
  //     TP_Allow_file.getline(u1, 256, '\n');
  //     TP_Allow_file >> Same_Particles_Flag;
  //     TP_Allow_file.getline(u1, 256, '\n');
  //     for (int i = 0; i < 6; i++) {
  //       int min, max;
  //       TP_Allow_file >> min >> max >> E_sign[i];
  //       // cout << "\t TP allow : " << min << "\t" << max << endl;
  //       // cout << A << "\t** " << endl;
  //       TP_allow_vec[0][i] = min;
  //       TP_allow_vec[1][i] = max;
  //       TP_Allow_file.getline(u1, 256, '\n');
  //     }
  //     TP_Allow_file.close();
  //     TP_allow_gather.push_back(TP_allow_vec);
  //     E_cut_gather.push_back(E_cut);
  //     E_sign_gather.push_back(E_sign);
  //     Same_Particles_gather.push_back(Same_Particles_Flag);
  //     single_t.print(TP_allow_vec);
  //   }
  // }
  cout << "Finished read data" << endl;
  // boost::mpi::broadcast(MainCom, TP_allow_gather, 0);
  // boost::mpi::broadcast(MainCom, E_cut_gather, 0);
  // boost::mpi::broadcast(MainCom, E_sign_gather, 0);
  // boost::mpi::broadcast(MainCom, Same_Particles_gather, 0);

  // triple_t.build_T(single_t,10);
  // triple_t.build_T(single_t,3,5,3,5,0,2);

  if (serial_flag_SC_cal == 1) {
    cout << "Start build_T" << endl;

    if (get_info.tp_info.Type == "abc") {
      int a_min = get_info.tp_info.abc_vec[0];
      int a_max = get_info.tp_info.abc_vec[1];
      int b_min = get_info.tp_info.abc_vec[2];
      int b_max = get_info.tp_info.abc_vec[3];
      int c_min = get_info.tp_info.abc_vec[4];
      int c_max = get_info.tp_info.abc_vec[5];
      triple_t.build_T(single_t, a_min, a_max, b_min, b_max, c_min, c_max);
    } else if (get_info.tp_info.Type == "E_max") {
      int Emax_t = get_info.tp_info.E_max;
      triple_t.build_T(single_t, Emax_t);
    } else {
      cerr << " Wrong @ build TP " << endl;
      cerr << "\t get_info.tp_info.Type : " << get_info.tp_info.Type << endl;
      exit(0);
    }
    // cout << "Finished build_T" << endl;
    triple_t.Select_FI_T(single_t, TP_allow_gather);
    // cout << "Finished Select_FI_T" << endl;
    if (MainCom.rank() == 0)
      triple_t.print("serial_dat/Output/TP_before_delete.dat");
    triple_t.Delet_Sort_T();
    // cout << "Finished Delet_Sort_T" << endl;
    triple_t.Delet_conf_T(single_t, E_sign_gather, E_cut_gather,
                          Same_Particles_gather);
    // cout << "Finished Delet_conf_T" << endl;
    triple_t.build_conf_T(single_t, E_sign_gather, E_cut_gather,
                          Same_Particles_gather);
    // cout << "Finished build_conf_T" << endl;
    triple_t.find_E_Jac_conf(E_Jac_conf_t);
    // cout << "Finished E_Jac_conf" << endl;
    if (MainCom.rank() == 0) {
      std::ofstream E_conf_file; 
      E_conf_file.open("data/Input/E_conf.Input");
      E_conf_file << E_Jac_conf_t.size() << "\t E_conf Num" << endl;
      for (int i = 0; i < E_Jac_conf_t.size(); i++) {
        E_conf_file << E_Jac_conf_t[i][0] << "\t " << E_Jac_conf_t[i][1]
                    << "\t E_F \t E_I" << endl;
      }
      E_conf_file.close();
    }
    if (MainCom.rank() == 0)
      triple_t.print_TPconf_info("serial_dat/Output/TPconf.dat",
                                 "serial_dat/Output/TPconf_E.dat",
                                 "serial_dat/Output/TPconf_E_tot.dat");
    triple_t.sort();
    // triple_t.Delet_conf_E_T();

    if (MainCom.rank() == 0)
      triple_t.print("serial_dat/Output/TP.dat");
    if (MainCom.rank() == 0)
      triple_t.print_channel("serial_dat/Output/TP_channel.dat");
    cout << "Start build_Tz" << endl;
    if (get_info.tp_info.Type == "abc") {
      int a_min = 2 * get_info.tp_info.abc_vec[0];
      int a_max = 2 * get_info.tp_info.abc_vec[1] + 1;
      int b_min = 2 * get_info.tp_info.abc_vec[2];
      int b_max = 2 * get_info.tp_info.abc_vec[3] + 1;
      int c_min = 2 * get_info.tp_info.abc_vec[4];
      int c_max = 2 * get_info.tp_info.abc_vec[5] + 1;
      triple_t.build_Tz(single_t, a_min, a_max, b_min, b_max, c_min, c_max);
    } else if (get_info.tp_info.Type == "E_max") {
      int Emax_t = get_info.tp_info.E_max;
      triple_t.build_Tz(single_t, Emax_t);
    } else {
      cerr << " Wrong @ build TP_Tz " << endl;
      cerr << "\t get_info.tp_info.Type : " << get_info.tp_info.Type << endl;
      exit(0);
    }
    triple_t.Select_FI_Tz(TP_allow_gather);
    triple_t.Delet_Sort_Tz();
    triple_t.Delet_conf_Tz(single_t, E_sign_gather, E_cut_gather,
                           Same_Particles_gather);
    triple_t.build_conf_Tz(single_t, E_sign_gather, E_cut_gather,
                           Same_Particles_gather);
    if (MainCom.rank() == 0)
      triple_t.print_TPconf_info_Tz("serial_dat/Output/TPconf_Tz.dat",
                                    "serial_dat/Output/TPconf_E_Tz.dat");
    // triple_t.Delet_conf_E_Tz();

    // triple_t.build_Tz(single_t,10);
    if (MainCom.rank() == 0)
      triple_t.print_Tz("data/TP_Tz.dat");
    if (MainCom.rank() == 0)
      triple_t.print_channel_Tz("data/TP_channel_Tz.dat");

    // triple_t.build_TzM(single_t,10);
    // triple_t.print_TzM("data/TP_TzM.dat");
    if (MainCom.rank() == 0 && serial_flag_SC_save == 1)
      ser_t.save(triple_t, "serial_dat/Output/TP_serial.dat");
  } else {
    ser_t.load(triple_t, "serial_dat/Output/TP_serial.dat");
    if (MainCom.rank() == 0)
      cout << "Finished Load triple_t" << endl;
  }

  Time_Count TC;

  TC.start();

  // cout << "Begin to build Jacobi basis" << endl;

  Jacobi jaco_b;
  Antisym antisym_t(sixJ_x, ninJ_x, hobra_t, jaco_b);
  Jacobi_NL jaco_nl;
  if (serial_flag_SC_cal == 1) {
    // jaco_b.build(3,2,3); //nmax l12max l3max
    // jaco_b.build_JT(2,2,3,3);
    // jaco_b.build_JT(8,3,3);

    jaco_b.build_JT(triple_t.Emax, triple_t.J2_max_TP + 2 * triple_t.Emax);
    if (MainCom.rank() == 0)
      jaco_b.print("serial_dat/Output/Jaco_unselect.dat");
    jaco_b.select_FI_T(triple_t.TPconf_E_tot);
    jaco_b.Delet_Sort_T(triple_t);

    jaco_b.sort();
    if (MainCom.rank() == 0)
      jaco_b.print("serial_dat/Output/Jaco_unselect2.dat");
    if (MainCom.rank() == 0)
      jaco_b.print_channel("serial_dat/Output/Jaco_channel_unselect2.dat");
    // jaco_b.print();

    double Anti_test =
        (1 - 2 * antisym_t.mat_cal(jaco_b.state[1][0], jaco_b.state[1][0])) / 3;
    if (MainCom.rank() == 0)
      cout << " \t ===+++=== TEST === Anti_test = " << Anti_test << endl;

    antisym_t.cal();
    if (MainCom.rank() == 0)
      antisym_t.print("serial_dat/Output/anti_vec.dat");
    if (MainCom.rank() == 0)
      jaco_b.print_channel("serial_dat/Output/Jaco_channel.dat");
    jaco_b.select();
    jaco_b.print("serial_dat/Output/Jaco_unselect3.dat");
    // jaco_b.print();
    // jaco_b.build_conf_T(triple_t.TPconf_E_tot);

    // exit(0);

    if (MainCom.rank() == 0) {
      cout << "Jacobi State Finished & Antisymmetry Time cost : " << TC.sec()
           << "s" << endl;
      cout << "Start to build Jaco_NL" << endl;
    }

    // jaco_b.print();

    // jaco_nl.build_NL(jaco_b, 7);
    if (MainCom.rank() == 0) {
      cout << "\t Triple Emax = " << triple_t.Emax << endl;
      cout << "\t Triple Emax_Final = " << triple_t.Emax_Final;
      cout << "\t Triple Emax_Init = " << triple_t.Emax_Init << endl;
    }

    jaco_nl.build_NL(jaco_b, triple_t);
    if (MainCom.rank() == 0)
      cout << "Finished build_NL" << endl;

    jaco_nl.select_FI_T(triple_t);
    // jaco_nl.Delet_Sort_T();
    if (MainCom.rank() == 0)
      cout << "Finished select_FI_T" << endl;

    vector<int> Jaco_used_Flag;
    jaco_nl.Delete_unused_Jaco_channel(jaco_b, Jaco_used_Flag);

    jaco_b.Delete_unused_Jaco_channel(Jaco_used_Flag);
    jaco_nl.sort();
    jaco_b.build_conf_T(triple_t.TPconf_E_tot);
    if (MainCom.rank() == 0)
      jaco_b.print_Jacconf("serial_dat/Output/Jacconf.dat");
    //    jaco_nl.print_short("data/Jaco_NL.dat");
    if (MainCom.rank() == 0)
      jaco_nl.print("serial_dat/Output/Jaco_NL.dat");
    if (MainCom.rank() == 0)
      jaco_nl.print_channel("serial_dat/Output/Jaco_NL_channel.dat");

    if (serial_flag_SC_save == 1 && MainCom.rank() == 0) {
      ser_t.save(antisym_t, "serial_dat/Output/anti_serial.dat");
      ser_t.save(jaco_b, "serial_dat/Output/jaco_serial.dat");
      ser_t.save(jaco_nl, "serial_dat/Output/jacoNL_serial.dat");
    }
  } else {
    ser_t.load(jaco_nl, "serial_dat/Output/jacoNL_serial.dat");
    ser_t.load(antisym_t, "serial_dat/Output/anti_serial.dat");
    ser_t.load(jaco_b, "serial_dat/Output/jaco_serial.dat");
    if (MainCom.rank() == 0)
      cout << "Finished Load antisym_t, jaco_b, jaco_nl" << endl;
  }

  if (MainCom.rank() == 0)
    jaco_b.print("serial_dat/Output/Jaco.dat");
  // jaco_b.check();
  if (MainCom.rank() == 0)
    jaco_b.print_channel("serial_dat/Output/Jaco_channel.dat");

  TC.start();

  if (MainCom.rank() == 0) {
    cout << "Jaco l12max = " << jaco_b.l12max << "\t l3max = " << jaco_b.l3max
         << "\t lmax = " << jaco_b.lmax << endl;
    cout << "Jaco n12max = " << jaco_b.n12max << "\t n3max = " << jaco_b.n3max
         << "\t nmax = " << jaco_b.nmax << endl;
    cout << "Jaco Emax_F = " << jaco_b.Emax_F
         << "\t Jaco Emax_I = " << jaco_b.Emax_I << endl;
  }

  if (State_Coeff_Break == 1) {
    all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
    cout << "Finished State & Coeff and Exit successfully on rank : "
         << MainCom.rank() << endl;
    exit(0);
    return 1;
  }
  // double hbar_omg;
  // vector<int> wave_const(13);
  // vector<vector<int>> E_Jac_conf;
  // E_Jac_conf = triple_t.E_Jac_conf;
  // if (MainCom.rank() == 0) {
  //   std::ifstream WaveInit_file;
  //   WaveInit_file.open("data/Input/Wave_init.Input");
  //   char u1[256]; // unuse
  //   for (int i = 0; i < 13; i++) {
  //
  //     int A = -1;
  //
  //     if (i == 1) {
  //       WaveInit_file >> hbar_omg;
  //     } else {
  //       WaveInit_file >> A;
  //     }
  //     // cout << A << "\t** " << endl;
  //     wave_const[i] = A;
  //     WaveInit_file.getline(u1, 256, '\n');
  //   }
  //   WaveInit_file.close();
  //   // WaveInit_file.open("data/Input/E_conf.Input");
  //   // int E_conf_num = 0;
  //   // WaveInit_file >> E_conf_num;
  //   // WaveInit_file.getline(u1, 256, '\n');
  //   // E_Jac_conf.resize(E_conf_num);
  //   // for (int i = 0; i < E_conf_num; i++) {
  //   //   E_Jac_conf[i].resize(2);
  //   //   int Ep, E;
  //   //   WaveInit_file >> Ep >> E;
  //   //   // cout << A << "\t** " << endl;
  //   //   E_Jac_conf[i][0] = Ep;
  //   //   E_Jac_conf[i][1] = E;
  //   //   WaveInit_file.getline(u1, 256, '\n');
  //   //}
  //   // WaveInit_file.close();
  //   string f_src = "data/Input/Wave_init.Input";
  //   string f_dst = "serial_dat/Output/Wave_init.Input";
  //   File_M.copy(f_src, f_dst);
  //   f_src = "data/Input/E_conf.Input";
  //   f_dst = "serial_dat/Output/E_conf.Input";
  //   File_M.copy(f_src, f_dst);
  // }
  // broadcast(MainCom, hbar_omg, 0);
  // broadcast(MainCom, wave_const, 0);
  // broadcast(MainCom, E_Jac_conf, 0);
  // const int lamb_cut = wave_const[0];
  // // int hbar_omg = wave_const[1]; // from int -> double
  // int reg_pow = wave_const[2];
  // int k_mesh = wave_const[3];
  // int k_scale = wave_const[4];
  // int p_mesh = wave_const[5];
  // int Emax_F = wave_const[6];
  // int Emax_I = wave_const[7];
  // int n_max = wave_const[8];
  // int l_max = wave_const[9];
  // int lbarS = wave_const[10];
  // int lbarS_G = wave_const[11];
  // int lbarS_P = wave_const[12];

  vector<vector<int>> E_Jac_conf;
  E_Jac_conf = triple_t.E_Jac_conf;
  const int lamb_cut = get_info.wave_info.Lambda_Cut;
  double hbar_omg = get_info.wave_info.hbar_omega;
  int reg_pow = get_info.wave_info.regulator_power;
  int k_mesh = get_info.wave_info.k_mesh;
  int k_scale = get_info.wave_info.k_scale;
  int p_mesh = get_info.wave_info.p_mesh;
  int Emax_F = get_info.wave_info.EF_max;
  int Emax_I = get_info.wave_info.EI_max;
  int n_max = get_info.wave_info.n_max;
  int l_max = get_info.wave_info.l_max;
  int lbarS = get_info.wave_info.lbarS;
  int lbarS_G = get_info.wave_info.lbarS_G;
  int lbarS_P = get_info.wave_info.lbarS_P;

  int EJC_num = E_Jac_conf.size();
  int EJC_max = 0;
  for (int i = 0; i < EJC_num; i++) {
    if (E_Jac_conf[i][0] > EJC_max)
      EJC_max = E_Jac_conf[i][0];
    if (E_Jac_conf[i][1] > EJC_max)
      EJC_max = E_Jac_conf[i][1];
  }
  string wave_serial_file =
      "../Wave_mpi_E/serial_dat/Wave_serial" + std::to_string(lamb_cut) +
      "_ho" + std::to_string(hbar_omg) + "_reg" + std::to_string(reg_pow) +
      "_k" + std::to_string(k_mesh) + "_ks" + std::to_string(k_scale) + "_p" +
      std::to_string(p_mesh) + "_EF" + std::to_string(Emax_F) + "_EI" +
      std::to_string(Emax_I) + "_EJC" + std::to_string(EJC_num) + "-" +
      std::to_string(EJC_max) + "_n" + std::to_string(n_max) + "_l" +
      std::to_string(l_max) + "_g" + std::to_string(lbarS) + "_G" +
      std::to_string(lbarS_G) + "_P" + std::to_string(lbarS_P) + "_sep.dat";

  // cout << "\t \t" << env_main.processor_name() << "\t " << wave_serial_file
  //     << endl;

  // if (MainCom.rank() == 0)
  cout << wave_serial_file << endl;

  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);

  // int test_press = 0;
  if (MainCom.rank() == 0) {
    cout << "Try & Test size? 1:yes" << endl;
    // cin >> test_press;
    // test_press = 1;
    cout << "Test press : " << Test_Press << endl;
  }
  boost::mpi::broadcast(MainCom, Test_Press, 0);
  if (Test_Press == 1) {
    if (MainCom.rank() != 0) {
      exit(0);
      return 1;
    }

    int size_trans = 0;
    /*for (int i = 0; i < trans_t.Tval_vec.size(); i++) {
      for (int j = 0; j < trans_t.Tval_vec[i].size(); j++) {
        for (int k = 0; k < trans_t.Tval_vec[i][j].size(); k++) {
          size_trans += trans_t.Tval_vec[i][j][k].size();
        }
      }
    }/**/

    cout << "TP all_size = " << triple_t.all_size_T << endl;
    cout << "Jacobi state all_size = " << jaco_b.all_size;
    cout << "\t Jacobi state Emax = " << jaco_b.Emax << endl;
    cout << "\t E_tot : " << endl;
    for (int k = 0; k < triple_t.TPconf_E_tot.size(); k++) {
      cout << "File No : " << k << endl;
      for (int i = 0; i < triple_t.TPconf_E_tot[k].size(); i++) {
        cout << "\t" << triple_t.TPconf_E_tot[k][i].E_F;
        cout << "\t" << triple_t.TPconf_E_tot[k][i].E_I << endl;
      }
    }
    cout << "----- E_cut_conf ------ " << endl;
    vector<vector<int>> E_cut_conf;
    triple_t.find_E_Jac_conf(E_cut_conf);
    cout << E_cut_conf.size() << "\t E_conf Num" << endl;
    for (int i = 0; i < E_cut_conf.size(); i++) {
      cout << E_cut_conf[i][0] << "\t " << E_cut_conf[i][1] << "\t E_F \t E_I"
           << endl;
    }
    cout << "Jacobi_NL state all_size = " << jaco_nl.all_size << endl;

    cout << " Trans-coeff size = " << size_trans << "\t with Memory ~ ";
    cout << size_trans * 8.0 / 1024.0 / 1024.0 << " Mb" << endl;

    cout << "nmax = " << jaco_b.nmax << "\t lmax = " << jaco_b.lmax
         << "\t J12_2max = " << jaco_b.J12_2max
         << "\t j3_2max = " << jaco_b.j3_2max << endl;
    Tpe_term tpe_test(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
    cout << "tpe_test start to init" << endl;
    tpe_test.init(jaco_b.Jaco_l123_info, jaco_b.Jaco_E_info, jaco_b.J12_2max,
                  jaco_b.j3_2max);
    tpe_test.get_Llbar_info();
    cout << "Ic1 size = " << tpe_test.Ic1_size << endl;
    cout << "Ic3 size = " << tpe_test.Ic3_size << endl;
    //<< "\t Ic3 lbarS = " << tpe_test.Ic3_lbar_max << endl;
    cout << "Ic4 size = " << tpe_test.Ic4_size
         << "\t Ic4 lbarS = " << tpe_test.Ic4_lbar_max << endl;
    Info_Wave IW_t;
    IW_t.n_max = jaco_b.nmax;
    IW_t.l_max = jaco_b.lmax;
    IW_t.lbar123_max = tpe_test.lbar123_max;
    IW_t.E_Jac_conf = E_Jac_conf;
    IW_t.Llbar_Ic_Reduced = tpe_test.Llbar_Ic_Reduced;
    IW_t.Llbar_Ic_Red_vec = tpe_test.Llbar_Ic_Red_vec;
    IW_t.l_vec = tpe_test.l_vec;
    IW_t.lbar_vec = tpe_test.lbar_vec;
    IW_t.l_state.clear();
    int size_l = tpe_test.l_state.size();
    for (int l = 0; l < size_l; l++) {
      l_struct_IW l_t;
      l_t.l12p = tpe_test.l_state[l].l12p;
      l_t.l3p = tpe_test.l_state[l].l3p;
      l_t.l12 = tpe_test.l_state[l].l12;
      l_t.l3 = tpe_test.l_state[l].l3;
      IW_t.l_state.push_back(l_t);
    }
    int size_lbar = tpe_test.lbar_state.size();
    IW_t.lbar_state.clear();
    for (int i = 0; i < size_lbar; i++) {
      lbar_struct_IW lb_t;
      lb_t.lbar = tpe_test.lbar_state[i].lbar;
      lb_t.lbar12 = tpe_test.lbar_state[i].lbar12;
      lb_t.lbar3 = tpe_test.lbar_state[i].lbar3;
      IW_t.lbar_state.push_back(lb_t);
    }
    Serial_SL_IW SL_IW;
    string serial_file_iw = "serial_dat/Output/info_wave.dat";
    SL_IW.save(IW_t, serial_file_iw);
    cout << "Finished Press Test" << endl;
    exit(0);
    return 1;
  }

  // ser_t.load(wave, wave_serial_file);
  Tpe_term tpe_t(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
  cout << "End to load wave Main on rank : " << MainCom.rank() << endl;
  if (serial_flag_WT_cal == 1) {
    TC.start();
    if (MainCom.rank() == 0) {
      cout << "\t ---- ---- Start Part : Wave & Tpe ---- ---- " << endl;
      cout << "Wave Fun Initial Time cost : " << TC.sec() << "s" << endl;
      cout << " Wave lbar_glbarS_max = " << wave.lbar_glbarS_max << endl;
      cout << " Wave lbar_PFglbarS_max = " << wave.lbar_PFglbarS_max << endl;
    }
    int mpi_size = MainCom.size();
    int mpi_rank = MainCom.rank();

    // all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
    // exit(0);
    /*wave.PFKglbarS_load_rank(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale,
                             p_mesh, Emax_F, Emax_I, n_max, l_max, lbarS,
                             lbarS_G, lbarS_P, mpi_rank, mpi_size, E_Jac_conf);
    cout << " rank = " << MainCom.rank() << "\t end to load\t "
         << env_main.processor_name() << endl;
    all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);/**/

    TC.start();

    int tpe_init_flag = 1;
    // exit(0);
    // if (tpe_init_flag == 1) {
    //   if (MainCom.rank() == 0) {
    //     cout << "lmax = " << jaco_b.lmax << "\t J12_2max = " <<
    //     jaco_b.J12_2max
    //          << "\t j3_2max = " << jaco_b.j3_2max << endl;
    //   }
    // } else {
    //   ser_t.load(tpe_t, "serial_dat/Output/Tpe_serial.dat");
    // }
    // tpe_t.build(jaco_b);
    // tpe_t.build_nlpnl(jaco_b);
    all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
    // cout<<"Test point 1"<<endl;
    tpe_t.init_II(jaco_b.Jaco_l123_info, jaco_b.Jaco_E_info, E_Jac_conf,
                  jaco_b.J12_2max, jaco_b.j3_2max);

    tpe_t.TPE_Build_MPI(E_Jac_conf);
    // tpe_t.Jlpl_build(tpe_t.l_state,tpe_t.l_vec,tpe_t.Jlpl);
    MainCom.barrier();
    all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
    cout << "Test point 1" << endl;
    // exit(0);
    // tpe_t.Jnlp_build(E_Jac_conf);

    // all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);

    // cout << " rank = " << MainCom.rank() << "\t start to build_MPI" << endl;
    // tpe_t.build_MPI();
    // cout << " rank = " << MainCom.rank() << "\t END to build_MPI" << endl;

    if (MainCom.rank() == 0) {
      tpe_t.print("serial_dat/Output/tpe.dat");
    }
    // tpe_t.print("data/tpe.dat");
    // tpe_t.C4_term(jaco_b.state[0][0], jaco_b.state[0][1]);
    if (MainCom.rank() == 0 && serial_flag_WT_save == 1)
      ser_t.save(tpe_t, "serial_dat/Output/Tpe_serial.dat");
  } else {
    ser_t.load(tpe_t, "serial_dat/Output/Tpe_serial.dat");
    cout << "Finished Load tpe_t" << endl;
  }
  if (MainCom.rank() == 0) {
    cout << "TPE term Time cost : " << TC.sec() << "s" << endl;
  }
  if (Wave_Tpe_Break == 1) {
    all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
    cout << "Finished Wave & Tpe and Exit successfully" << endl;
    exit(0);
    return 1;
  }

  //-----------------------------Trans ----------------------------------//

  if (MainCom.rank() == 0) {
    cout << "Start Init TransCoeff" << endl;
    TransCoeff trans_t(sixJ_x, ninJ_x, hobra_t, jaco_nl, single_t, triple_t);
    // jaco_nl.print(2);
    // triple_t.print(2);
    // trans_t.build_JTE(1,1,2);
    // double test=trans_t.cal(jaco_nl.state_NL[0][0],triple_t.state[0][0]);
    // cout<<"T test = "<<test<<endl;

    if (serial_flag_OCP_cal == 1) {
      trans_t.build();
      if (serial_flag_OCP_save == 1) {
        ser_t.save(trans_t, "serial_dat/Output/trans_serial.dat");
      }
      trans_t.print("serial_dat/Output/trans_vec.dat");
    } else {
      ser_t.load(trans_t, "serial_dat/Output/trans_serial.dat");
    }
    cout << "T-coefficients Time cost : " << TC.sec() << "s" << endl;

    //---------Cont term, Ope term & Package ---------------//
    wave.init_single(jaco_b.nmax, jaco_b.lmax);
    Cont_term cont_x(LECs, wave, sixJ_x, jaco_b);
    if (serial_flag_OCP_cal == 1) {
      // cout << " nmax : " << jaco_b.nmax << endl;
      cont_x.init(jaco_b.nmax); //( n_max )
      // cout << "good 111 ?" << endl;
      cont_x.build();
      // cout << "good 2  ?" << endl;
      cont_x.print("serial_dat/Output/cont_term.dat");
      if (serial_flag_OCP_save == 1) {
        ser_t.save(cont_x, "serial_dat/Output/Cont_serial.dat");
      }
    } else {
      ser_t.load(cont_x, "serial_dat/Output/Cont_serial.dat");
      cont_x.print("serial_dat/Output/cont_term.dat");
    }
    cout << "Contact term Time cost : " << TC.sec() << "s" << endl;
    TC.start();
    // jaco_b.print(5);

    Ope_term ope_x(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
    if (serial_flag_OCP_cal == 1) {
      ope_x.init();
      // ope_x.cal(jaco_b.state[5][68],jaco_b.state[5][51]);
      ope_x.build();

      ope_x.print("serial_dat/Output/ope_term.dat");
      // ope_x.print();
      if (serial_flag_OCP_save == 1) {
        ser_t.save(ope_x, "serial_dat/Output/Ope_serial.dat");
      }
    } else {
      ser_t.load(ope_x, "serial_dat/Output/Ope_serial.dat");
      ope_x.print("serial_dat/Output/ope_term.dat");
    }
    cout << "OPE term Time cost : " << TC.sec() << "s" << endl;

    Package pack(hbar_omg, triple_t, jaco_nl, jaco_b, trans_t, antisym_t,
                 cont_x, ope_x, tpe_t);

    cout << "Begin to build_JacAvec : " << endl;
    if (serial_flag_OCP_cal == 1) {
      pack.build_JacAvec();
      // pack.print_JacAvec("data/JacA_vec.dat");
      pack.build();
      pack.print("serial_dat/Output/package.dat");
    }
    cout << "Start to save package serial" << endl;
    if (serial_flag_OCP_save == 1) {
      ser_t.save(pack, "serial_dat/Output/Package_serial.dat");
    }
    // pack.print_nnlosat(single_t, "data/sat/Navr_ord_states.dat",
    //                   "data/sat/package_sat.dat");
    cout << "Package Time cost : " << TC.sec() << "s" << endl;
    cout << "Total time cost : " << TC_main.sec() << "s" << endl;
    int print_flag = 0;
    if (serial_flag_OCP_save == 1 && print_flag == 1) {
      string filename_TP_mat_ALL = "serial_dat/Output/package_ME_ALL.dat";
      std::ifstream NO_file;
      NO_file.open("data/Input/Model_space.Input");
      if (!NO_file) {
        std::cout << "File does not exist : "
                  << "data/Input/Model_space.Input" << std::endl;
        exit(0);
      }
      vector<vector<int>> model_vec(5);
      for (int i = 0; i < 5; i++) {
        char u1[256]; // unuse
        int A, B;
        NO_file >> A >> B;
        model_vec[i].push_back(A);
        model_vec[i].push_back(B);
        NO_file.getline(u1, 256, '\n');
      }
      NO_file.close();
      int sp_core = model_vec[0][0];
      int sp_core_Tz = model_vec[0][1];
      int sp_P_beg = model_vec[1][0];
      int sp_P_beg_Tz = model_vec[1][1];
      int sp_P_end = model_vec[2][0];
      int sp_P_end_Tz = model_vec[2][1];
      int sp_Q_beg = model_vec[3][0];
      int sp_Q_beg_Tz = model_vec[3][1];
      int sp_Q_end = model_vec[4][0];
      int sp_Q_end_Tz = model_vec[4][1];
      std::ifstream OC_file;
      std::string file_oc = "data/Input/ReOrder_SP.Input";
      OC_file.open(file_oc);
      if (!OC_file) {
        cout << "File does not exist : " << file_oc << endl;
        exit(0);
      }
      int oc_num;
      char u_c[256];
      OC_file >> oc_num;
      OC_file.getline(u_c, 256, '\n');
      vector<int> order_change(oc_num);
      for (int i = 0; i < oc_num; i++) {
        int x, y;
        OC_file >> x >> y;
        cout << x << "\t " << y << endl;
        OC_file.getline(u_c, 256, '\n');
        order_change[x - 1] = y - 1;
      }
      OC_file.close();
      pack.print_TP_ALL(filename_TP_mat_ALL, single_t, sp_core_Tz, sp_P_beg_Tz,
                        sp_P_end_Tz, hbar_omg, reg_pow, lamb_cut, order_change);
    }
    cout << " ==== ALL FINISHED ==== " << endl;
    // TC.start();
  } /**/
}
